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ABSTRACT 

We have extended the spectral dynamics formalism introduced by Binney & Spergel, 
and have implemented a semi-analytic method to represent regular orbits in any poten- 
tial, making full use of their regularity. We use the spectral analysis code of Carpintero 
& Aguilar to determine the nature of an orbit (irregular, regular, resonant, periodic) 
from a short-time numerical integration. If the orbit is regular, we approximate it by 
a truncated Fourier time series of a few tens of terms per coordinate. Switching to a 
description in action-angle variables, this corresponds to a reconstruction of the under- 
lying invariant torus. We then relate the uniform distribution of a regular orbit on its 
torus to the non-uniform distribution in the space of observables by a simple Jacobian 
transformation between the two sets of coordinates. This allows us to compute, in a 
cell-independent way, all the physical quantities needed in the study of the orbit, in- 
cluding the density and in the line-of-sight velocity distribution, with much increased 
accuracy. The resulting flexibility in the determination of the orbital properties, and 
the drastic reduction of storage space for the orbit library, provide a significant im- 
provement in the practical application of Schwarzschild's orbit superposition method 
for constructing galaxy models. We test and apply our method to two-dimensional 
orbits in elongated discs, and to the meridional motion in axisymmetric potentials, 
and show that for a given accuracy, the spectral dynamics formalism requires an order 
of magnitude fewer computations than the more traditional approaches. 

Key words: galaxies: kinematics and dynamics - celestial mechanics, stellar dynam- 
ics 



1 INTRODUCTION 

In the construction of dynamical models for galaxies by 
Schwarzschild's (1979) method, one tries to match the den- 
sity distribution and the photometric and kinematic ob- 
servations with weighted contributions of individual orbits 
(e.g., Rix et al. 1997; Cretton et al. 1999). Therefore, the 
method requires the precise knowledge of the intrinsic prop- 
erties of the orbits, such as the density and the line-of-sight 
velocity distributions (losvd's). The traditional method for 
computing the orbital properties is to integrate the equa- 
tions of motion by numerical means. Then the spatial den- 
sity distribution of the orbit follows from calculating the 
fraction of time spent in each cell of a grid in configura- 
tion space after a sufficiently long integration period. Other 
properties follow by considering similar grids in phase-space. 
This method can be visualized as simply marching along 
the orbit, and dropping balls in a grid of buckets at regular 
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time intervals. It is conceptually similar to a Monte-Carlo 
integration, and can be applied to any orbit. However, this 
approach does not take advantage of the regular behavior of 
many orbits in the potentials relevant for galaxies: it simply 
treats the regular orbital motion as a structureless collection 
of independent points in phase-space. 

The motion in a regular orbit is quasi-periodic, and 
can be expressed as a Fourier series expansion in action- 
angle variables. This series expansion represents the under- 
lying invariant orbital torus (e.g., Arnold 1989), and can be 
reconstructed from a normal orbit integration (Binney & 
Spergel 1982, 1984), or from generating functions (McGill 
& Binney 1990, Binney & Kumar 1993, Kaasalainen & Bin- 
ney 1994, Kaasalainen 1994). It is also possible to integrate 
the orbit directly in the angle variables (Ratcliff, Chang & 
Schwarzschild 1984), although published applications have 
been restricted to two-dimensional motion. Laskar (1993) 
introduced the so-called frequency map analysis to study 
the nature of orbits. This has been used to study orbits in 
galactic potentials (Papaphilippou & Laskar 1996, 1998; Val- 
luri & Merritt 1998). Carpintero & Aguilar (1998, hereafter 
CA98) developed a fully automated technique to classify or- 
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bits based on the commensurability of the peak frequencies 
in the Fourier spectra of the orbit. 

In this paper we extend this previous work as follows. 
From a standard numerical integration of the orbit over a 
few hundred orbital periods, and after a detailed spectral 
analysis based on the method of CA98, we obtain a semi- 
analytic expression for the orbit. This can be considered 
as fitting the underlying orbital torus to the results of a 
short numerical integration, and so obtaining an approxi- 
mate expression for the orbit valid for all time (i.e., fully 
phase-mixed). It is then straightforward to derive all the 
required orbital properties, including the density in config- 
uration space and the losvd's, with high precision, and we 
describe in detail how to do this. While there are many pa- 
pers on orbital tori reconstruction in the literature, we are 
not aware of any published formalism for projecting the tori 
to observable space to get the LOSVD of an orbit. Here we 
implement our method on two-dimensional integrable and 
non-integrable potentials, and compare the results with tra- 
ditional straight numerical integrations. We also show how 
to generalize the results to axisymmetric three-dimensional 
potential. 

This paper is organized as follows. After reviewing some 
basic results of spectral dynamics in we present the ex- 
tended formalism in §[| In ^ we test the formalism on a 
two-dimensional separable (Stackel) potential, and then ap- 
ply it to Binney's logarithmic potential. We describe the 
application to the motion in the meridional plane of an ax- 
isymmetric potential in §|B| and summarize our conclusions 
in §| 



2 INTRODUCTION TO SPECTRAL 
DYNAMICS 

For the convenience of the general reader, we here summa- 
rize some results on the structure of phase-space, and intro- 
duce the terminology and notations of action-angle descrip- 
tion of orbits. 



2.1 Phase-space structure, action-angle variables 
and base frequencies 

In an n-dimensional potential (cases of interest are n — 2 or 
3), an orbit can be characterized according to its m integrals. 
In a time-independent potential, the energy is a conserved 
quantity and is always an integral for an orbit, leading to 
m > 1. When the orbital motion in phase-space (of dimen- 
sion 2n) is constrained by the conservation of m integrals of 
motion, the orbital manifold on which the orbit evolves has 
dimension 2n — m. The orbit is said to be regular if m > n, 
and then travels on a manifold of dimension < n; otherwise 
(m < n), the orbit is said to be irregular (or chaotic), and 
move in a region of phase-space of dimension > n (but still 
< 2n - 1). 

In the special case m — n, the orbital manifold of a 
bound orbit is topologically equivalent to an n-torus, the in- 
variant torus, embedded in phase-space (e.g., Arnold 1989; 
see Fig. |l|), and the orbital motion is quasi-periodic. Not be- 
ing further constrained by additional integrals of motion, the 
trajectory is dense (or ergodic) on the torus, and the orbit 
is said to be open. On the other hand, if m > n, the orbit is 




Figure 1. A regular open orbit spiraling on its invariant torus. 



further constrained, and will not cover densely its torus even 
after infinite time. This happens when there exist resonances 
between the fundamental frequencies of two or more angle 
variables describing the orbit (cf. following section). Thus, 
we speak of non-resonant (m = n) and resonant (m > n) 
orbital tori; only the former form a set of non-zero measure 
in phase-space. If the orbit is fully resonant, i.e., m = 2n— 1, 
it closes on itself only after a finite number of turns around 
its torus: such an orbit is periodic. 

A regular orbit is confined to its invariant n-torus, so 
we can look for canonical coordinates (Q, P) adapted to the 
torus geometry, such that the momenta P are constant on 
the torus, and their conjugate variables Q form a natural 
coordinate system for the region of phase-space occupied by 
the tori. Since these coordinates are canonical, Hamilton's 
equations state: 



Q 



dH 
dP' 



dH 
9Q' 



(1) 



where H is the Hamiltonian of the orbit. Since we want 
P to be constant, this implies that H — H(P), i.e., that 
H is cyclic in Q, and thus that dH/dP = cte = a;. We 
can therefore integrate the remaining equation of motion 
straightforwardly: Q(t) = w t + Q(0). 

Choosing an appropriate scaling such that Q is 2-7T- 
periodic, we obtain the so-called action-angle variables: the 
angles ip = Q are natural coordinates for the invariant torus, 
and allow one to distinguish individual points on it; the ac- 
tions J = P label this torus with respect to the tori of other 
orbits. They can be seen as related to the different radii of 
the tori. 

If the frequencies (actually pulsations) uJi (i = 1, . . . , n) 
are mutually incommensurable, they are said to be the 
72bf = n base frequencies of the orbit; the incommensu- 
rability ensures that the orbit covers densely its torus and 
therefore is open. If there exist a linear commensurability be- 
tween two or more uii, then the orbit is further constrained 
on the torus, and the number of incommensurable base fre- 
quencies is reduced to 71bf < n. In the extreme case, if the 
orbit is closed, only one base frequency remains. In the other 
extreme case, the base frequencies are not well-defined for ir- 
regular orbits, since these are not constrained to a n-torus; 
however, according to the practical definition of the base 
frequencies, they will be characterized by hbf > n (CA98). 
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2.2 Description of an orbit in angle variables 

We now turn to the two-dimensional case, n = 2, and restrict 
ourselves to regular open orbits, which leads to bbf = n — 2. 
The other regular orbits (with «bf < 2) can be considered 
as a degenerate subset of this case, in the sense that the 
resonances lead to invariant tori of zero measure in phase 
space. 

A regular open orbit in a two-dimensional potential is 
characterized by its two base frequencies ui = 2ttv. The 
angle variables <p £ [0, 2ir[ are then given by 



cp(t) = wi, 



(2) 



where we have chosen the zero of time such that <p(0) = 0. 
We can then write the position r = (x, y) and the velocity 
v = f at time t as 



r(t) = r(<p(t)) 



and 



>(t) = v(<p(t)). 



(3) 



The quasi-periodicity of the orbit allows us to write the or- 
bital motion in a Fourier form (Binney & Spergel, 1982): 



x(t) = E^(!,m) cos ( w (l,m)*+X(!,m)), 

l.m 

V(t) = J2 Y (l,m) COS ( UJ (Lm,)t + i'(l,, n )), 



(4) 



where (I, m) C 1? are pairs of integers, X(i,m) an d ipa im ) 
are constant phases, and the frequencies u)a tm \ axe linear 
combinations of the base frequencies: 



W( ij>n ) = luii + mu>2 = A 



(5) 



and we have used the shorthand A = (l,m). Alternatively 
we can write 



x(t) — ^2 x >> cos(A • w t + xa), 



(6) 



general to a still-to-be-defined number of solutions ipg 
such that: 



r(<Po ) = r o, 



i = 1, 



,N. 



(10) 



The total number N of solutions is finite (see, e.g., Fig. ^). 
We will refer to the N different single-valued functions in- 
volved in the definition of <p by ^ , such that r(<p^) can 
be inverted. 

The non- uniqueness of the solution can be understood, 
in the case of <p(r) for instance, by the fact that an orbit 
usually has a finite number of ways to reach a given point 
ro, either from a different part of the trajectory, or along 
the same part of trajectory but in different directions in a 
case of an orbit without a given sense of propagation (such 
as a box orbit). Accordingly, for a given orbit we can pre- 
dict the general number of solutions iV (excepting the de- 
generate cases) of eq. ([lo]). Specifically, N = 2 for an open 
two-dimensional tube orbit, since the orbit can usually ac- 
cess a given point along two paths with a given direction of 
propagation imposed by the fixed sense of rotation. N — 1 
on the boundaries of the orbit. For a flat box orbit each of 
the two paths through a given point can be traveled in both 
directions since there is no fixed sense of rotation, so that 
N — 4. On the outer boundaries of the box, degeneracy re- 
sults in A^ = 2, or even iV = 1 in the corners of the orbit, 
which are accessible in only one way (the one perpendicular 
to the zero- velocity curve; see Ollongren 1962). Each time 
a resonance brings two distinct parts of the orbit on top of 
each other, it is easy to see that the number of solutions in 
the zone of overlap will be doubled. For instance, we expect 
N = 8 in the zone of overlap of an open 3 : 2 ('fish') orbit 
(Miralda-Escude & Schwarzschild 1989; CA98). 



for the position and 

v x (t) = — ^2 A ■ uj X\ sin(A ■ wi + x\) 



(7) 



for the velocity, and similar expressions for y(t) and v y (t). 

Switching from time t to angle variables up by means 
of eq. (^|), we obtain the torus expression of the motion (in 
what we will refer to as the (^-space hereafter), in the sense 
that angle variables are natural coordinates for the invariant 
torus on which the regular orbit evolves: 

A 

Vxjtp) = - / ] X\ X • w sin(A • (p + x\), (9) 

A 

and similar expressions for y(<p) and v y (tp). 
3 EXTENDED FORMALISM 



3.1 



(^-inversion 



The first step for calculating the density p at a point ro = 
(xo,yo) in configuration space is to find the corresponding 
actions and angles at that point. However, while r(tp) and 
r(<p) (cf. eqs (g) and (|^)) are single-valued functions, the 
inverse functions <p(r) and tp(r) are generally multi- valued. 
For instance, a given position ro on the orbit corresponds in 
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3.2 Orbital density 

In configuration space, we can define the orbital density p 
at point ro by writing the mass of the orbit enclosed in the 
element of volume dr = dxdy around ro as 



dm = p(ro) dr. 



(11) 



Equivalently, in y>-space, we can define the density p^ at 
a point (po such that the mass enclosed in the element of 
volume d<f> = dtpi dipi around tpo is: 



dm' = p P (<po) dip. 



(12) 



If we choose to normalize the total mass of the orbit to 1, 
then J dm = J dm' = 1. 

An open orbit, i.e., an orbit with two incommensurable 
base frequencies, will eventually fill its torus uniformly, so 
that the density p^ in <^>-space is constant ('time averages 
theorem', BT87). Mass normalization then gives: 



Pv = 



Jd<t> (2tt) 2 



(13) 



If the orbit is closed, i.e., the orbit has only one base fre- 
quency, all the other frequencies in eq. (^) being commen- 
surable with it, then the time averages theorem does not 
apply, since the orbit is confined to a closed spiral on its 
torus by a new integral of motion. 

Using the single- valued functions giving the solu- 
tions ip^ of eq. ([w]) , we can relate the element of volume dr 
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around the point ro to the elements of volume d(f>^ around 



(0 



d0 (l) = 



dr 



dr = 



dr 



J, 



(0 : 



where we have defined the Jacobian Jr as in eq. (|1 



T (i) 



Or 



dip 



(14) 



(15) 



We can then link the density p to the constant density p v 
given by eq. ( |l3| ) by properly relating eqs ([y]) and jl^): 



dm = dm' 1 ' 



(16) 



with dm' 2 ' = p v {jp^ ) d<f>^ , since the mass lying in dr 
around a given point is the sum of the masses lying in the 
various elements d(f>^ contributing to dr. This gives: 



dm = 



dr 



1 



(17) 



which, by comparison with eq. (|ll[), leads to the final ex- 
pression: 

1 - 1 

the Jacobians being computed at the N solutions ip^ 
of eq. (0). The various J^ 1 ' do not have to be the same, 
since the functions (p^' are different. 



3.3 Boundary of the orbit 

The boundary B of an orbit can be considered as the lo- 
cation of the points leading to the degeneracy of solutions 
of eq. ([lo|). At a point rg of B, we can say that if ip& is a 
solution of r(ifi) = rg, then r (ipg + Sip) is also a solution to 
the second-order in ||<fys||; we have: 



r(ip B + Sip) 

if and only if 
dr 



r{(fis) + 



dr 



Jv(<pb) 



dp 



dp 

= r B + 0(\\5ip\ 
0. 



x Sip 



IMI 



2\ 



(19) 



We then reach another equivalent definition of the boundary 
B, as the location of the points where at least one Jacobian 
J r vanishes (and possibly more, since there is a degeneracy 
of the solutions at this point). 

At the boundary, the orbit is divergent in the density 
(cf. eq. (|is|)), but this divergence is integrable since the total 
mass of the orbit is finite. 



3.4 Orbital distribution function 

The above considerations, related to equations analogous 
to (^) and ([]), directly lead to an expression of the or- 
bital distribution function (df) /(ro,i>o), by considering 



^(2rr) 2 J^'j as the contribution to the local density p from 
solution i among N: 

1 N 1 

f(ro,v ) = £ JL x S (v - «,(„«)) , (20) 



where <p 



(0 



, N are the solutions to eq. (110). Given 



the relation (|18[), this expression is consistent with the basic 
relation: 



p( r ) = J f{r,v) 



dv, 



(21) 



where the velocity is an indirect function of the position. 
Eq. ( |20[ ) gives the df at any position ro. 

There are several variations for expressing the phase 
space density, depending whether we are primarily inter- 
ested in the distribution of the position or of the velocity. 
For generalization, let us consider a vector k from a half- 
subspace (of dimension n = 2) of phase-space (of dimen- 
sion 2n = 4), and i from the complementary half-subspace. 
The previous situation corresponded to k — r = (x, y) and 
l = v = (v x , v y ), but we could have, as we shall see in §3.5, 
k = (:c s ky, vios) and l = (xios, Wsky)- We now look for an 
expression of the df f(Ko,to) from k. 

In the same way as in the special case k = r described 
one can show that: 



in 63.2 



11 



(2tt) 2 ^ tW 

v ' i=l J n 



where this time, we have defined Jr as 



J, 



(i) 



dip 



(22) 



(23) 



and with ip^ being the multiple solutions of the equation: 
i=l,...,M. (24) 



k(Vo ) = K o, 



The number M of solutions does not have to be the same as 
the number N of solutions to eq. (|l^) . 

3.5 Line-of-sight velocity distribution 

Eq. (|2o| ) provides the df f(r, v) at any point of phase-space 
accessible to the orbit, and hence allows computation of all 
the orbital dynamical quantities, directly related to the df, 
in particular the line-of-sight velocity distribution (losvd). 

Once we have chosen a line-of-sight (los), associated 
with coordinates (x',y') (y' along the LOS and x' being in 
the plane of the sky), we can express the LOSVD V~P x i{v y i) 
(the velocity profile) by the usual formula: 

VP.'(v)= / dy' fdv x ,f(r',v'). (25) 

J los J 

We see that, for this purpose, we need an expression for 
the df as a function of £ = (x',v y /). According to eq. (pj|), 
we can write: 



VP 



1 M 1 
= (2^)2 53 JW' 



(26) 



where the Jacobians are to be computed at the M so- 
lutions of the equation: 
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C(v>&°) = Co 



(27) 



We can also write down the surface brightness /u(x') as: 
(i(x) = / dy / dv'f(r', v') 

JIOB J 

= / p(r')dy'= [vP x ,(v)dv, 

JIOB J 

which gives two equivalent ways to compute these quantities, 
either from the density, or from the LOSVD. 

3.6 Actions 

The formalism we use here, based on an expression of the 
dynamical quantities of the orbit in terms of the angle vari- 
ables (p, is well-suited for the computation of the actions. 
An open two-dimensional regular orbit is characterized by 
exactly two integrals of motion, which can be for instance 
the two actions J\ and J%, that can be defined by (BT87): 



v ■ dr, 



1,2. 



(28) 



(p 4 e[0,2ir[ 



Using relations (Q) and (B), we obtain after some algebra: 



J, 



G) 



(29) 



which is very similar to eq. (16) in Binney & Spergel (1984). 
As explained by these authors, one is free to define the 'real' 
(useful) actions J r and J a as any linear combination of the 
previous quantities, so that it matches any natural require- 
ment such as J r = for closed long-axis orbits and J a = 
for closed loop orbits (BT87) . In our case, with the CA98 ex- 
traction procedure, we found that the linear transformation 
to apply was (cf. Fig. [Io| : 



Jr = 
Ja — 



and 



Jr 
J a 



Ji 

J-2 



2J 2 
Ji - 



for box orbits, 



J> 



for loop orbits. 



(30) 



(31) 



This choice guarantees a continous action space (cf. Binney 
& Spergel 1984; de Zeeuw 1985). 



4 2D NUMERICAL IMPLEMENTATION 

We have numerically implemented the formalism described 
above, and in §kj| describe tests on orbits in a two- 
dimensional separable potential for which all relevant quan- 
tities are known analytically. We th en consider a non- 
integrable logarithmic potential in §4.2. 



4.1 Sridhar &: Touma potential 

We first study an orbit in a two-dimensional separable 
potential, described by Sridhar & Touma (1997; see Ap- 
pendix M). This potential is of Stackel form in parabolic 
coordinates, and allows a complete analytic study of the 




Figure 2. Fourier spectra (amplitude vs. frequency) for the test 
orbit. The selected frequencies for the decomposition (squares) 
are those of amplitude Ai larger than F max[Ai] in each coordi- 
nate, indicated by the long-dashed lines, with F = 10 — 3 (giving 
15 peaks in x and 26 in y). The two base frequencies are marked 
by stars. 



orbits (in particular computation of the boundaries and ex- 
pression of the surface density), which are all regular, and 
usually denoted as banana orbits. 

The illustrations in this section are made for a slightly 
cusped potential (a = 0.5, with notations of Appendix^), 
and a test-orbit with energy E — 2.39 and second integral 
I2 = 0.38, integrated on 4096 points over 123 periods with 
a 7/8 th order Runge-Kutta integrator (Fig. 0). 



4-1.1 Spectral analysis and time- dependent expressions 

The Fourier spectra x{y) and y(v) of the time series of the 
coordinates x(t) and y(t) of a quasiperiodic orbit, obtained 
in the usual way by numerical integration, consist of discrete 
lines whose frequencies are linear combinations of the base 
frequencies. We modified an algorithm provided by CA98 
for the extraction of the base frequencies, in order to give 
all the quantities defining the Fourier series (Q). The CA98 
method computes the Fourier transform of the coordinates 
of an integrated orbit, identifies the peaks in the Fourier 
spectra and extracts the corresponding frequencies. It looks 
for the base frequencies, and linearly decomposes the se- 
lected frequencies over these BFs (see Fig. ^|). The output 
then consists of: the base frequencies w, and for each peak 
frequency in each coordinate, the linear decomposition over 
the base frequencies (e.g., A = (I, m)), and the associated 
amplitude (X\) and phase (xa)- 

The CA98 algorithm provides a finite number of terms 
in the quasi-periodic expansion (^J). The selection is made 
by keeping, in each coordinate, only frequencies with ampli- 
tude greater than a fraction F of the greatest one (typically 
F — 10 -2 -10 -3 , depending on the required accuracy), giv- 
ing a finite number of terms in each direction (typically be- 
tween 10 and 20). The numerical approximations obtained 
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-1 
x 



Figure 3. Comparison between time expression from numerical 
integration (left) and spectral expansion (right) of an orbit in the 
Shridhar & Touma potential. The spectral time-dependent expan- 
sion (^) was obtained with F = 10~ 3 (resulting in 15 frequencies 
in x and 26 in y) at the same time-steps as the integrated orbit. 
The short-dashed line is the Zero- Velocity Curve (zvc), and the 
long- dashed curves indicate the boundaries of the orbit calculated 
by using the separable nature of the potential (cf. Appendix 
The orbit was launched from the ZVC (square). The boundaries of 
the reconstructed orbit are barely distinguishable from the exact 
boundaries, indicating the level of precision of the expansion. 



by the truncation of the series coincide with the actual ex- 
pansion (ft) up to the required accuracy (Fig. 0). 




,t/2 * 3tt/2 2h 

Pi 



Figure 4. Contours of x((p) (solid lines) and y(ifi) (dotted lines), 
for the test orbit reconstructed with F = 10 — 3 . The levels x = 
xo = 0.5 and y = yo = (thick lines) intersect at the TV = 4 values 
of (fi(ro) (dots). Successive initial conditions for the root-finding 
routine (squares) lead to different solutions (straight lines), until 
it finds the expected number of solutions. 



4-1.2 ip- expressions 

The time-dependent expansion (g) provides a semi-analytic 
expression for a regular orbit. It can be useful for storage 
purposes, or for extrapolating orbital motion to times much 
longer than the actual integration time (e.g., in order to 
reduce the numerical noise of a grid-based density computa- 
tion). However, time is not a good variable for a regular orbit 
that will eventually fill up completely its invariant torus, and 
the angle variables <p should be used instead. This leads to 
the expressions (^) and (^) for x(tp) and y(tp), respectively. 

T he p roblem is then to recover the N different solutions 
of eq. (lid) for a given position ro. First, we have to find the 
expected number N of solutions. This is inferred from the or- 



bit classification, according to the criteria described in §3.1 



We then use standard root-finding routines (in our case, the 
NAG routine c05pbf) from different initial conditions until 
we have found the JV solutions, or until a maximum number 
of trials, indicating a degenerate case with fewer solutions 
(see Fig. Q). Interestingly, there is point-symmetry in the so- 
lutions. For example, for the Sridhar & Touma banana orbit 
in Fig. ^, the center of the symmetry is at ip = (w, 7r). In 
fact this is true for any box orbit in a general non-rotating 
potential as well. This is because we can always choose the 
corner of a box orbit as the initial conditions, i.e., v — (p — 
at t — 0. Then for any box described by r(t), we can gener- 
ate the new orbit described by r(—t) by reversing the arrow 
of the time, and the new orbit simply repeats the old one 
with r(— t) — r(t): reversing has no effect on the orbit be- 
cause the initial velocity is zero. This time-symmetry means 
ip [2tt] and —tp [2tv] describe the same point in the configu- 
ration space, where the modulus brings the phase angle to 
the default interval [0, 2ir]. Thus the solutions come in pairs 
at <p and 2n — (p with point-symmetry around (tv, 7t). The 



above arguments do not apply to loop orbits because there 
is no point in a loop where the velocity is zeroQ. 



4-1.3 Density maps and velocity profiles 

Eq. ([l8]) allows the computation of the orbital density at 
any point ro once we have determined the N solutions of 
eq. (llOj) by the appropriate root-finding routine. Since the 
test orbit was integrated in a Stackel potential, the derived 
density distribution can be compared to the analytic expres- 
sion available for this case (cf. Appendix |a|) . Fig. ^ shows 
contours of the two densities side by side with the same lev- 
els. The spectral density was obtained from an expansion 
truncated in amplitude with F = 10 -3 , giving 15 terms in x 
and 26 in y. The indicated masses, which should be unity by 
construction in both cases, are computed from a numerical 
integration over the density grid; it can be therefore used as 
a rough estimate of the precision of the density computation. 

As expected, the more terms are kept in the series (1), 
the more accurate the representation of the orbit (Fig. 5). 
With F = 10 -2 , the orbit is represented by an expansion 
with 4 terms in x and 12 terms in y, and the density is 
recovered to ~ 20%. When F is set to 10~ 3 (Fig. |), the 
expansion contains 15 terms in x and 26 in y, and the density 
is accurate to ~ 3% in the inner parts of the orbit, and 
degrades very close to the boundaries of the orbit due to the 
divergence of the density there. This precision in density 
can be compared with the usual com putation of the density, 
directly during integration (see §4.3). 



t Reversing the arrow of time from a initial point at say, x = 
0, y = 1 of a loop, will lead to a counter-loop with identical shape 
indeed, but for which ip and 2tx — ip correspond to two different 
points in configuration space, with (x(—t), y(—t)) = (—x(t), y(t)). 
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Figure 5. Upper panel: Orbital density distribution in the Shrid- 
har &; Touma potential computed from spectral theory (right) 
compared with the exact expression (left). The spectral den- 
sity ( |li| ) was obtained with F = 10~ 3 on the same grid as the 
Stackel density (step 5 gr id = 10 — 2 ). The density levels range from 

to 2 in 11 steps. The total mass of the orbit is normalized to 

1 by definition in each case. Lower panel: Ratio between spectral 
density and Stackel density, in the case frac = 10 — 3 . The con- 
tour levels range linearly from 0.9 to 1.1 in steps of 0.01 (thick 
line = 1). The increased discrepancy close to the boundary of the 
orbit is due to the divergence of the density and its associated 
numerical imprecision. 



Stackel 



Spectral 




frac = lCr B 
M. - 0,973 



Figure 6. Position-velocity diagram (pvd), i.e., density in the 
plane (x,y), for the test-orbit, computed from Stackel expression 
(left) and from spectral theory with F = 10 -3 (right). The con- 
tours range linearly from to 0.4 in 9 steps. Since the orbit is 
a banana orbit, hence a boxlet orbit without a definite sense of 
rotation, the PVD is symmetric with respect to y = 0. 



Using eq. (|26|) , we can also compute the velocity profile 
along any line-of-sight, and gather in a position- velocity di- 
agram (pvd) all the velocity profiles corresponding to paral- 
lel LOS. This is the projected density of the phase-space orbit 
on the plane (x',v y i) (where x' and y' are the coordinates 
running perpendicular and parallel to the LOS), in the same 
way as the 'density' corresponds to the projected density of 
the phase-space orbit on the plane (x,y). A cut of this pvd 



along a given x' would give the LOSVD along a line-of-sight 
running parallel to the y'-axis and going through x' . Fig. ^ 
presents side by side the PVD of the test-orbit for a LOS 
parallel to the y-axis, computed from the Stackel formalism 
(left) and from the spectral theory (right), using the same 
spectral expansion as befor e. Pre cision of the computation 
is as expected (~ 10%, see §4.2.1). 



4.2 Logarithmic potential 

We now consider a non-integrable potential of astronomical 
interest, the logarithmic potential (BT87, p. 126): 



$(x,y) = i«g In 



Rl + x 2 + K 



(32) 



with Wo the circular speed at large radius, R c the core radius, 
and q the axial ratio of the equipotentials. Positivity of the 
associated density requires l/v2 < q < 1. This potential 
admits two major families of orbits (BT87), the box orbits 
and the flat-tube or loop orbits, accessible from different 
initial conditions. For the illustrations, we consider the case 
vo = R c = 1 and q — 0.9. 



4-2.1 Density maps and velocity profiles 

Two orbits of the same energy E — 1.629 were integrated: a 
box orbit over 128 periods, and a loop orbit over 167 periods 
(cf. Fig. Q). This energy gives a zvc such that xzvc /Rc = 5, 
where xzvc corresponds to the intersection of the zvc with 
the x-axis. 

From eq. (p"s|), we then compute the density map of 
each orbit from a spectral analysis set by F — 10 -3 , giving 
position expansions with 11 terms in x and 17 in y for the 
box orbit, and 13 terms in x and 12 in y for the loop orbit 
(right panels of Fig. ^). 

Fig. |^ presents the pvd of the two previous orbits for 
a LOS parallel to the y-axis. The 'wiggles' that can be seen 
in the middle of the loop orbit pvd (Fig. ^ lower panel) are a 
numerical artifact, at the ~ 10% level. The origin of this can 
easily be understood: for an orbit in logarithmic potential 
([32]), the characteristic length scale, velocity and frequency 
are of the order xzvc, «o and |w| ~ vo/xzvc- When we 
truncate the expansion (^) of the coordinate x(<p) according 
to the parameter F, we make a relative error at the level 
X\/xzvc < F with the coordinate, where X\ is the am- 
plitude of the truncated term. Meanwhile, the correspond- 
ing term in the velocity expansion (^) has an amplitude 
A ■ u> X\ , where the frequency of the expansion term A • w 
increases linearly as we go to higher and higher order terms. 
So the corresponding error with the velocity is at the level 
A ■ us X\/vq ~ |A|Xa/:tzvc < |A|.F. So while the same trun- 
cation reproduces the orbital position and position-related 
quantities such as the density within a few percents, it may 
significantly over- or under-estimate the high frequency vari- 
ations in the orbital velocity and for velocity-related quanti- 
ties such as losvd's. A way to correct for this effect would be 
to truncate the different expansions according to the higher 
amplitude X\ or A a; X\ at a given order (Papaphilippou & 
Laskar 1996). While truncating according to the amplitudes 
X\ is straightforward, truncating according to A ■ u> X\ re- 
quires a priori knowledge of the BF decomposition (A,w), 
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Figure 7. Upper panel: Integrated box orbit (left) and spectral 
density (right) computed from an expansion truncated at F = 
10 — 3 on a grid of step <5 gr id = 10 — 2 . The contours range from 
to 0.1 in 11 steps. The short- dashed line is the ZVC, and the 
long-dashed line corresponds to the core of the potential. Lower 
panel: Same as previous, but for the y > 0-side of the loop orbit. 
The contours range from to 5.10 - 2 in 11 steps. 



E = 1 629, v = 1. R p = 1, 0.90 




12 3 4 5 

x (dy/dt>0) 



Figure 9. Surface of section (a;, x), y > computed from 60 reg- 
ular open orbits of energy E = 1.629, for a truncation fixed by 
F = 10 -3 . Only the low -resonance orbits which do not overlap 
over themselves, have been kept. The empty areas correspond 
either to stochastic areas or 'high'-resonance islands. Dots corre- 
spond to the construction of the surface of section of each orbit 
during its short numerical integration (see text). 




Figure 8. Upper panel: PVD for the box orbit, with contours 
ranging linearly from to 0.3 in 11 steps. Notice the symmetry 
of the PVD with respect to y = 0. Lower panel: Same as previous, 
but for the loop orbit. The contours range linearly from to 0.2 
in 11 steps. Since the loop orbit has a definite sense of rotation 
(here, anti-clockwise), the PVD only has reflection symmetry with 
respect to the point (x,y) = (0,0). 



which generally would mean a double pass, and hence a sig- 
nificantly increased number of computations. We have not 
followed this route here, as this additional accuracy will play 
only a minor role in most practical applications, such as the 
computation of the moments of the LOSVD. 



4-2.2 Surface of section and action space 

We have integrated a complete library of orbits in the loga- 
rithmic potential (|32|). All the orbits have the same energy 
E — 1.629, and are integrated over typically ~ 150 peri- 
ods with 4096 time-steps. The 140 initial conditions (ICs) 
are distributed according to: 40 ICs evenly spaced in angle 
on the zvc, and thus with zero initial velocity so that they 
are boxes or irregular orbits, and 100 ICs evenly spaced in 
radius along the short axis (x = and y > 0), with an ini- 
tial velocity vector perpendicular to the axis (v x > and 
v y — 0). Those with y < yb, where yt is the amplitude of 
the last stable t/-axis oscillation, are boxes or irregular orbits. 
The remainder are loops. With the exception of the closed 
loop at the chosen energy, these all cross the i/-axis perpen- 
dicularly in two points. We keep only those that cross at 
apocenter, and furthermore remove all high-order resonant 
orbits, leaving a total of 60 regular, open, 'low-resonance' 
boxes and loops. The spectral analysis of each integrated 
orbit is carried out with F = 10 -3 , giving typically ~ 13 
terms in the expansions along x and y. 

We constructed the surface of section (x,x), y > for 
this library by solving for the velocity (x,y) at each point 
(x,y = 0) of an orbit. The result is shown in Fig. ^. Since 
we kept only low-resonance orbits, this surface of section 
does not show any resonant island. However, the discarded 
orbits leave empty areas in the surface of section between 
the regular loops and boxes. For comparison, we also show 
in in Fig. [] the surface of section from the short orbital in- 
tegration integration (~ 150 periods). As already noted in 
§4.2.1, the velocity is less accurately retrieved in the spec- 
tral process than the position, which results in some small 
differences (few percents) in the surface of section. 

Using expression fl29|), we computed the actions as- 
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Figure 10. Action space J a , vs. J r computed from the same or- 
bital set as Fig. |jj and with the same truncation of the expansion. 
Note the natural distinction occurring between box orbits (solid 
squares) and the loop orbits (open squares). The two straight se- 
quences meet in the middle as in Fig. 3-25 of BT87. The gaps in 
the surface of section plot are seen here as well. 



Figure 11. Comparison, for the x > 0, y > part of the box 
orbit of Fig. [j| between the spectral density (upper left panel and 
dotted line in other panels) and the 'bucket' densities (solid line in 
other panels) for different integration times, related to the number 
of points. In both cases, contour levels range logarithmically from 
1.35 10~ 4 to 4 10 — 4 in 4 steps, and the step of the grid used for 
the computation of the density is 5 SI ^ = 0.1. The mean accuracy 
p. is defined in the text. 



sociated with an orbit, and constructed the action-space 
(JJ Tl Ja) of the mono-energetic library (Fig. [loj) . This vi- 
sualization of the library in action-space is of particular in- 
terest: since the action-angle variables are canonical, and 
their integration is particularly simple (§2.1), equal volume 



in action space is directly associated to equal volume in 
phase-space (Binney & Spergel 1984). This provides a two- 
dimensional representation of the four-dimensional phase- 
space. As expected, a natural dichotomy between box orbits 
(solid squares) and loop orbits (open squares) appears in 
this space. 



4.3 Numerical discussion 

The traditional way of computing the orbital density, just 
as any other dynamical quantity such as the LOSVD, is the 
following: along the orbital motion (during the numerical 
integration), and for a predetermined grid in configuration 
space, one counts up the number of integration points falling 
in a given cell (e.g., Rix et al. 1997). This gives the total 
fraction of time spent in this cell, which is proportional to 
the local density after a long enough integration time. 

For a sufficiently fine grid, this process of 'dropping balls 
in buckets' can be considered as a random process, and the 
statistical noise associated with a count of points in a cell 
is ~ \/iV, leading to a precision in density of order l/x/TV. 
The orbit in Fig. ^| was integrated for 4096 points, and the 
density represented on a grid with step <5 gr id — 10" 2 . As- 
suming that the integration points are approximately uni- 
formly distributed over the whole area accessed by the orbit, 
roughly 2 x 0.5 = 1 or l/(5g rid = 10 4 cells, the average num- 
ber of points per cell is ~ 0.5, leading to a very low precision 
in density (> 100%). A precision of the same order of the 
one obtained with spectral dynamics (^ 5%), would require 
integrating ~ 800 times longer, in order to get ~ 400 points 
per cell. 

Fig. |ll] and [li] compare in a more quantitative way the 
spectral density (dotted line) and the 'bucket' densities (solid 
line) computed for different integration times for the box and 




4 1 



Figure 12. Same as Fig. [ll| but for loop orbit of Fig. [?]. Contour 
levels range logarithmically from 1.1 10 — 4 to 2.3 10 -4 in 4 steps. 



loop orbits shown in Fig. tn. The estimated mean accuracy jX 
is denned as fl = (l/\/Ni), where Ni is the number of balls 
in bucket i, and (•) is the average over all the non-empty 
cells. As expected, the mean accuracy is directly related to 
the number of integration points used during the computa- 
tion of the density. It can be clearly seen in both cases that 
4096 points along the orbit, typically needed for the torus 
reconstruction from spectral analysis, are highly insufficient 
for computing a proper 'bucket' density. 

Fig. |l^ compares the different mean accuracies obtained 
in terms of computing power. The CPU-time is expressed in 
arbitrary units. Each density was computed on the same 
computer, and the time used by each process was renor- 
malized to the shortest one (a few seconds on a standard 
Linux station) . The spectral mean accuracy is estimated to 
be ~ 3% (cf. §4.1.3). We conclude that the spectral density 
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Figure 13. Mean accuracy of the density computations vs. com- 
puting power needed for the computation, for the previous box 
(solid symbols) and loop (open symbo s) orbits. Each dot on the 
lines corresponds to a case in Fig. [ll or |l^, while the isolated 
symbols give an estimate of the spectral density mean accuracy 
(~ 3%). 



is almost one order of magnitude more accurate than the 
'bucket' density for the same computing power. 




Figure 14. Three-dimensional axisymmetric orbit integrated in 
the meridional plane (upper part) and its associated meridional 
density computed from spectral analysis (lower part). 



5 AXISYMMETRIC POTENTIALS 

Recent work on modeling the internal dynamics of elliptical 
galaxies, aimed at measuring masses of central black holes 
or the distribution of extended dark matter, compare stellar 
absorption line kinematics with fully general three-integral 
multi-component axisymmetric models. While in some cases 
semi-analytic approaches can be used (e.g., Matthias & Ger- 
hard 1999), many studies use a variant of Schwarzschild's or- 
bit superposition technique (e.g., van der Marel et al. 1998; 
Cretton & van den Bosch 1999; Gebhardt et al. 1999; Cret- 
ton, Rix & de Zeeuw 1999) . In this approach, orbits are cal- 
culated numerically in a chosen potential, and their prop- 
erties are projected onto the plane of the sky. These are 
then combined to reproduce the observed surface density 
and losvd's. The current implementations use the 'bucket' 
method for calculating the orbital properties. 

Motion in a three-dimensional axisymmetric potential 
can be reduced to a two-dimensional problem by exploit- 
ing the conservation of the z-component of the angular mo- 
mentum L z (e.g., BT87), and describing the motion in the 
meridional plane (R, z) (the angular variable (f> being cyclic). 
The properties of a regular orbit in such a potential can thus 
be computed with much improved accuracy by means of a 
minor adaptation of the spectral dynamics formalism de- 
scribed in the above, and we present the relevant equations 
here, together with some examples. 

The illustrations are made for a test-orbit with energy 
Eo — —0.8 and angular momentum L z — 0.2 integrated in 
the meridional plane of a three-dimensional axisymmetric 
coreless logarithmic potential: 



*(fl,z) = i«gln [R 2 + ^ 



(33) 



with q = 0.9 and vo = 1. Hereafter, we use a = i — n/2, 
where i is the usual inclination angle of a galaxy (i — for 
a face-on galaxy). 



meridional plane 
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Figure 15. Notations used in the text for the conversion of two- 
dimensional to three-dimensional densities. We use a = i — 7r/2, 
where i is the customary inclination of the galaxy (> when seen 
from above). 



5.1 Spatial density 

If the density of an orbit in the meridional plane is pm(R, z), 
the associated spatial density is simply 



p(R,z,(p) = 



Pm(R,z) 



(34) 



2-kR ' 

where the factor 2ir ensures the normalization, so that 
/ pmARAz = J pRdRdzdcf) = 1. 

Fig. [l4| shows an example of the density of an orbit in 
the meridional plane, reconstructed with the spectral for- 
malism. 



5.2 Projected intensity 

Cartesian coordinates (x' , y') on the plane of the sky re- 
late to spatial cylindrical coordinates (R, (j>, z) as follows (see 
Fig. El: 
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Figure 16. Sky density computed for the test-orbit seen edge- 
on (a = 0, upper panel) and slightly overhead (a = 10°, lower 
panel). Dotted line: ZVC in the meridional plane. Note how the 
divergence of the density at the boundaries of the orbit in the 
meridional plane translates in sky density Heavy line: geometrical 
location of the 'corners' of the orbit. The numerical line-of-sight 
integration precision is set to 10 — 2 . 



x' = R cos cf>, 

y = z cos a ~ R sin (f> sin a, 



(35) 



where a defines the orientation of the line of sight. Assuming 
a 7^ ±7r/2 (i.e., the galaxy is not seen face-on), we find 



R 



■x',v{<t>) 

Zx',y>{<t>) 



\x'\/ cos 4>, 

j/'/cosa + \x'\ tanevtanc 



(36) 



Furthermore, since dz' 
dz'= W, . 



-dy/cosa and dy = Rd(j>/ cost 



(37) 



We can then compute the intensity on the sky I a (x',y') = 
J los p(R, z, <f>) dz' for x' 7^ 0, to obtain: 



I<x{x',y) 



2tt 



71/2 pM{R X >,y>(<f>),Z x ,, y ,{cf>))< 

. w / 2 cos a cos 



(38) 



The integrand is not singular at the limits of integration, 
since an orbit reaches a maximum radius J? max in the merid- 
ional plane, corresponding to a maximum angle ci m ax = 
arccos(7?o/i?max) (for Ro < i? ma x; otherwise, we know that 
1 = 0). Equivalently, the orbit reaches a minimum radius 
i? m in, corresponding to the angle 4>min = arccos(J?o/^min) 
(or if Ro > -Rmin). We can therefore restrict the integration 
domain to the segments [— </> ma x, — <^min] U \ 
Similarly, for x' = 0, we obtain 



I a (x' = 0,2/') = 



dz' 



(39) 



where we have defined p., , (z) = pm(z' cos a, — — ±z' sin a). 

I 1 V ' ' r il + v i C os a > 

Fig. |lfj shows the upper-right quadrant of the projected 
density on the plane of the sky of the test orbit when seen 




Figure 17. Elements of geometry for the computation of the 
velocity along the LOS. 



edge-on (a = 0) or slightly overhead (a = 10°). As the 
computation of the projected intensity involves an extra nu- 
merical integration, the result is more sensitive to numeri- 
cal errors in the determination of the orbital density in the 
meridional plane. Fig. |l^ has been obtained with an inte- 
gration relative precision set to 10 -2 . 



5.3 Line-of-sight velocity distribution 

The velocity «i os along the LOS (cf. Fig. ^) is given by 
L 



^ios = I —- cos <p + vr sin < 



cos a + v z sin a, 



(40) 



where vr and v z are the velocity components in the merid- 
ional plane, and v$ = L z /R. The LOSVD can then be com- 
puted as follows: 



VP. 



'(Wlos) 



(277) 

where the Jacobian 

(i) ~ d{x',y',vi os ) 



1 M 1 

— T — 



J(i)- 



J 



is to be evaluated at the i-th solution ip^' = 
i = l,..., M, of the system of equations: 

= z(ip^) cos a — R(ip^) sin</>sin a, 



(41) 



(42) 



V\ os = 



L, 



R( V ^)) 



cos <j> + vr(ip^) sin 4>J cos a 

+ v z ((p^) sin a 



(43) 



The reason that we need only to solve for two action angles 



and 



a(0 



for the motion in the meridional plane is that 



(i) 

the third angle <f>. a is simply the azimuthal angle 
is related to the two action angles <p^ by 



= arccos 



R(<p 



(01 



which 



(44) 



Fig. [l8| shows the position-velocity diagram of the test- 
orbit as seen with a 'virtual' long-slit placed along the x'-axis 
(that is y' = 0). By symmetry, we can restrain the study to 
the part x > 0, where the velocities are expected to be 
mostly positive since L z > 0. 



6 CONCLUSIONS 

We have used the spectral analysis code provided by Carpin- 
tero & Aguilar (1997) to determine the properties of orbits 
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Figure 18. PVD for a 'virtual' long-slit placed along the x'-axis 
(y 1 = 0) of the test-orbit seen edge-on (a = 0). The intensity at 
a point (x 1 ,y' = 0,v\ OB ) is computed by means of eq. 

in elliptic disks and in axisymmetric potentials based on nu- 
merical integration over a modest number of orbital periods. 
We have approximated regular orbits by a truncated Fourier 
time series of a few tens of terms per coordinate, and have 
reconstructed the underlying invariant torus by computation 
of the associated action-angle variables. We have used this 
to relate the uniform distribution of a regular orbit on its 
torus to the non-uniform distribution in the familiar space 
of observables by a simple Jacobian transformation between 
the two sets of coordinates. This approach is by no means 
new, but we have extended the published formalism to in- 
clude calculation of the quantities used in the construction 
of dynamical models by Schwarzschild's method, in particu- 
lar the orbital density and the observed losvd's for elliptic 
disks and for axisymmetric models. 

In the standard implementation of Schwarzschild's 
method, orbital properties are recorded on a grid of cells 
in configuration space (x,y,z) and in the space of observ- 
ables (x' , y' , v\ OB ) (e.g., Cretton et al. 1999). The number of 
cells is generally taken to be of the order of 1000-2000, and 
orbital densities and kinematic properties are calculated by 
numerically integrating an orbit for a very long time, so that 
each cell is crossed at least 100 times. Averaging of the time 
spent in the cell, or of the velocity distribution in the cell, 
then results in a ~ 10% accuracy on the orbital density and 
the observed kinematics. By contrast, the spectral method 
allows the invariant orbital torus to be described by a few 
tens of terms, and allows computation of all the physical 
quantities in a cell-independent way with much increased 
accuracy. At a given accuracy, the spectral method requires 
less cpu time than the traditional approach, at least for the 
two-dimensional cases we have investigated. The resulting 
flexibility in the determination of the orbital properties, and 
the drastic reduction of storage space for the orbit library, 
provide further significant improvements in the practical ap- 
plication of Schwarzschild's method. 

The algorithm devised by CA98 uses only the orbital 
position for spectral analysis. It truncates the orbital Fourier 
series according to a criterion on the spectra of the coor- 
dinates, so that the resulting Fourier series approximation 
and torus reconstruction is less accurate in velocity space. 
Laskar's (1993) naff algorithm uses the complete phase- 
space position, and may therefore be more efficient in the 
determination of the base frequency of an orbit and in the 
torus reconstruction. However, the precision obtained with 
CA98 seems adequate for the practical goals of the method. 

The method employed here to compute orbital prop- 



erties works best for the main families of regular or- 
bits. The higher-order resonances (e.g., Miralda-Escude & 
Schwarzschild 1989) can be found by the CA98 algorithm, 
but the torus reconstruction becomes more difficult (e.g., 
Kaasalainen 1995). This suggests that the main advantage 
will lie in the construction of dynamical models that have 
a fairly regular phase space, such as models with central 
density profiles that are only moderately cusped. 

Irregular orbits cannot be represented by our Fourier 
series expansions. In principle it is possible to obtain the 
properties of these orbits by subtracting the contributions of 
the regular ones from an f(E) or f(E, L z ) component (e.g., 
Zhao 1996; Cretton et al. 1999), which can be represented as 
the sum of all regular and irregular orbits at that E (and L z ). 
Use of the smooth regular orbits and the f(E) and f(E, L z ) 
components as building blocks in Schwarzschild's method 
therefore means that the irregular orbits are included as 
well. 

Application of spectral dynamics to compute the ob- 
servable properties of regular orbits in genuinely three- 
dimensional potentials is possible in principle. Construction 
of triaxial dynamical models that fit observed the kinemat- 
ics of early-type galaxies may well benefit from the same 
approach, although the rich phase space, and the presence 
of many irregular orbits, may complicate its practical appli- 
cation (e.g., Hafner et al. 1999). 

It is a pleasure to thank Luis Aguilar and Daniel Carpin- 
tero for making available their spectral dynamics code, and 
to thank Wyn Evans for constructive comments. 
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APPENDIX A: SRIDHAR & TOUMA 
POTENTIAL 

The test-orbit used in §^[]] was integrated in the two- 
dimensional separable potential introduced by Sridhar & 
Touma (1997, ST). We collect here a number of properties 
of the orbits in this potential. 

The ST potential is defined by 



$(r,(9) = r a [(!• 



+ (1 - COS0) 



(l + c) 



(Al) 



It is of Stackel form in parabolic coordinates £ = r(cos8 + 
1) > and rj = r(cos 8 — 1) < 0, in which it takes the form 



$ = 



F+(0 F_(tj) 



with F+(0 = 2£ 1+Q and F.( V ) = -2|77| 1+Q . 
In addition to the energy 

2 



E = 



(?P5 " VPr,) + 



(A2) 



(A3) 



every orbit has a second isolating integral of motion I2 , given 
by 



J 2 = 2£pf-£E + F+(£) 



■r,E + F-(r,). 



(A4) 



As eq. (Al) shows, the ST potential is scale-free, so 



that the shape of an orbit is determined only by the value of 
the second integral, while its scale is related to the energy. 
Two orbits of opposite I2 are symmetric with respect to the 
:r-axis; we consider hereafter only the case I2 > 0. 

Defining the functions I+(rf) = — r\E + F-(rf) and 
J-(f) = — £E + F+(£) for a given energy E, the bound- 
aries of an orbit with second integral I2 (> 0) are given by 
the roots of the equations I+(rj) = I2 (2 roots r\\ and 772) 
and J-(£) = h (one root £0; see e.g., Fig. |). 

The orbital df of an orbit with energy Eo and second 
integral Jo can be written (as a consequence of the strong 
Jeans theorem) as: 



f Eo ,i (E, h) oc S(E - E ) S(h - Jo). 



(A5) 



The orbital surface density is as usually given by 
PE ,i (r,6) = J f Eo ,i d 2 v. We have d^drjdp^dp^ = drdv, 
which translates to: 

j2 o , _ 1 dj£yrf) djp^pr,) 



r d(r,0) d{E,I 2 ) 



dEdh, 



with 

d(r,8) 
and 



2r sin 8, 



(A6) 



(A7) 



djp^Prj) = 1 1 

9(E,h) 2 y /(i + (r,)-Io)(Io-I-(0)' 
Therefore, the unnormalized orbital density p* is 

Pk ^' V) = 2 V (/ +W _/ )(/ -/_ (e ))' 

while the total mass of this unnormalized orbit is found to 
be: 



(A8) 



(A9) 



A I 



B , -To 



PEoJo 1 



V2 -#d, r° 

vi v S v Jo V^s/Si 



+ 



1)2 



d.77 



vi VVv St, Jo 



to 



d£ 



(A10) 
(All) 

. (A12) 



with shorthand notations S v = ~ ^0 an d = Jo — 

J_(£). This double integration can be carried out easily by 
numerical means. 
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